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We develop a quantum algorithm to solve combinatorial optimization problems through quantum simulation 
of a classical annealing process. Our algorithm combines techniques from quantum walks, quantum phase 
estimation, and quantum Zeno effect. It can be viewed as a quantum analogue of the discrete-time Markov 
chain Monte Carlo implementation of classical simulated annealing. Our implementation requires order 1 / \f5 
operations to find an optimal solution with bounded error probability, where 8 is the minimum spectral gap of 
the stochastic matrix used in the classical simulation. The quantum algorithm outperforms the classical one, 
which requires order 1/5 operations. 
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I. INTRODUCTION 

Combinatorial optimization problems (COPs) such as the 
traveling salesman problem are important in almost every 
branch of science, from computer science to statistical physics 
and computational biology |1]. A COP consists of a family 
of instances of the problem; each instance is an optimization 
problem, to minimize (or maximize) some objective function 
over a finite set S of d elements, called the space of states. 
This space may have additional structure (e.g., it may be a 
graph), allowing the definition of a notion of locality; and 
the set of objective functions may have special properties de- 
pending on the particular COP. In general multiple local min- 
ima may be present. Finding a solution by exhaustive search 
is hard in general, due to the large size of the search space. 
Therefore, more efficient optimization approaches are desir- 
able. Efficiency is typically quantified in terms of how the 
resources needed to find the optimum scale with the instance 
size, which is typically polynomial in log d. 

Simulated Annealing (S A) is a possible generic strategy for 
solving a COP |2||. The idea of SA is to imitate the process 
undergone by a metal that is heated to a high temperature and 
then cooled slowly enough for thermal excitations to prevent 
it from getting stuck in local minima, so that it ends up in 
one of its lowest-energy states. In SA, the objective func- 
tion plays the role of energy, so the lowest energy state is 
the optimum. This process can be simulated using different 
techniques; we focus on discrete Markov chain Monte-Carlo 
(MCMC). These methods are often used to numerically obtain 
properties of, for example, classical physical lattice systems in 
equilibrium |3fl. The general idea of MCMC is to stochasti- 
cally generate a sequence of states via a process that converges 
to a target probability distribution. This is the Boltzmann dis- 
tribution at the low final temperature in the case of SA. The 
efficiency of the method relies on the fact that, in general, only 
a small proportion of states contribute significantly to the de- 
termination of properties in equilibrium. Therefore, if a good 
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state-generating rule is chosen, the MCMC algorithm can ex- 
plore the most relevant states only, outperforming exhaustive 
search. 

One way to estimate the implementation complexity of SA 
using MCMC is to count the number of times that the state- 
generating rule must be executed (i.e., the number of gener- 
ated states) in order that the desired distribution is reached 
within an acceptable error. This complexity, denoted by Ms a, 
is of order 0(log(d/e 2 )/<5) (see Sec. |H]). Here, S is the min- 
imum spectral gap of the stochastic matrices used to gener- 
ate states for the COP via MCMC 01, while e is the error 
probability, that is, the probability that the final state sampled 
via this process is not a solution (not in the set §o of optimal 
states). Ideally, Ms a is insignificant compared to the size of 
the state space. This is the situation, for example, when com- 
puting physical properties of the Ising spin model using the 
Metropolis rule yfl. In this example Ms a is known to be of 
order 0(N 2 ) for a system of N spins, while the state space 
dimension is d = 2 N . Nevertheless, Ms a can increase rapidly 
with N if the interaction strengths are made random J5D, mak- 
ing the problem intractable in general. In this case, this is 
due to the gap 5 becoming exponentially small in N (instance 
size). Therefore, finding new methods with better scaling in 
S, yielding speedups over SA, is of great importance. 

Quantum mechanics provides new resources with which to 
attack these optimization problems 10, 0. @l- Quantum com- 
puters (QCs) can theoretically solve some problems, includ- 
ing integer number factorization and search problems, more 
efficiently than today's conventional computers 19(]. Still, 
whether a QC can solve all COPs more efficiently than its 
classical counterpart is an open question. In this paper we 
show that QCs can also be used to speed up the simulation 
of classical annealing processes. That is, we present a new 
quantum algorithm that can be seen as the quantum analogue 
of SA using MCMC, but for which the number of times that 
the state-generating rule is called (Mqsa) is greatly reduced 
to 0(\og 3 {d/e 2 )/(VSe 2 )), to achieve error bounded by e, in 
a single run. This speed-up is most significant for hard in- 
stances where 6 -C 1. Our quantum simulated annealing 
algorithm (QSA) is constructed using ideas and techniques 
from quantum walks fiol [Till and quantum phase estima- 
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tion lfl2l [l3ll . Th e Q SA also exploits the so-called quantum 
Zeno effect [lj], in which after Q = 0(1/ At) mea- 
surements of a quantum system at short time-intervals At the 
state is collapsed onto the ground state with total probability 
1 - O(At). 

This paper is organized as follows. First, in Sec. [II] we de- 
scribe the implementation of SA using discrete-time MCMC, 
and in Appendix lAl we derive a rate at which the tempera- 
ture of a classical system can be lowered to assure conver- 
gence to the set of ground states. To do this we adapt the 
results obtained for the continuous-time case in Ref. H. The 
rate that we obtain is similar to the one in Ref. (To\ for those 
cases where S decreases exponentially with the problem size 
(cf. Ref. 01). In Sec. [Ill] we describe a quantization of a re- 
versible Markov chain in terms of quantum walks. Our quan- 
tization is a similarity-transformed version of the one used in 
Refs. iTToL rTvfl to speed up search problems. It constructs, from 
the transition matrix of the Markov chain, a unitary operator 
acting on a set of quantum states corresponding to the classi- 
cal ones. In Sec. [lV]we describe our QSA and obtain the cor- 
responding implementation complexity, exhibiting a quantum 
speed-up with respect to classical SA. Since our QSA makes 
calls to the phase estimation algorithm, we describe phase es- 
timation in Appendix 151 Finally, we present the conclusions 
in Sec|V] 



general the annealing schedule may strongly depend on (3k- 
In our case the overall implementation complexity of the al- 
gorithm with constant A/3 is of the same order as for a general 
annealing schedule, so the analysis below is valid for both sit- 
uations. 

We choose A/3 = 0(5/ Em), where 5 is the minimum 
spectral gap of the matrices M(/3k) at inverse temperature 
Pk = kA/3, and Em '■= max a \E[a]\. In Appendix lAl we 
show that for f3f = C(7 _1 log(d/e 2 )), the probability of not 
ending in a solution is no greater than e [see Eq. ( IA19M . 7 is 
the spectral gap of E. The implementation complexity of SA 
is then given by P = (3f /A/3. We obtain 

^ = 003^ = 0^^) (1) 

for a success probability greater than 1 — e. The dependence 
of Ms a on (5 _1 is characteristic of Markov processes and, al- 
though Eq. (T]i only gives an upper bound on the resources 
required for the implementation of SA, such a dependence on 
the spectral gap may be unavoidable lfl8ll . 

Remarkably, a similar algorithm implemented on a quan- 
tum computer has a reduced implementation complexity for 
those hard instances where i « 1. This is described in the 
following sections. 



H. SIMULATED ANNEALING AND MONTE-CARLO 
TECHNIQUES FOR MARKOV PROCESSES 

We consider the simulation of a classical annealing process 
via MCMC, and give annealing rates such that the final sam- 
pled state is almost certain to be in the set §0 of optimal solu- 
tions to a COR To do this, we first need a formulation of the 
COP in terms of an equivalent problem in which §0 consists of 
the states that minimize some real-valued cost function E on 
the state space. Usually, E is regarded as the energy function 
of a classical system S, so the optimal solutions to the COP 
are represented by the ground states of S. For concreteness, 
we sometimes think of S as defined on a lattice with N ver- 
tices, having a finite state space {a} of size d = 0(exp(N)). 

A ground state can be reached by annealing slowly enough, 
starting with S at sufficiently high temperature. The MCMC 
simulation of this process, described in terms of the inverse 
temperature (3 = 1/T, begins by sampling a state from 
the uniform distribution. The annealing process is determined 
by a choice of an annealing schedule, i.e. a finite increas- 
ing sequence /3i < 02 < ■ ■ ■ < (3p, and by a sequence of 
transition rules {M(j3k)}- Each M((3k) is a stochastic matrix 
whose elements m aa i ((3k) are transition probabilities from a 
to a'. M((3t) is chosen to have the Boltzmann distribution at 
(3k as its unique equilibrium distribution. 

At each step k, a new state is stochastically generated 
from er^" -1 ) according to the transition probabilities M (/?&). 
The annealing schedule is chosen to give an acceptable upper 
bound e on the probability of error (of not ending up in So). 
For simplicity, we consider an annealing schedule such that 
A/3 = Pk - (3k-i < 1 is constant, and thus (3 f = PA/3. In 



III. QUANTUM WALKS AND ERGODIC MARKOV 
CHAINS 

Discrete-time quantum walks were introduced as the quan- 
tum analogues of classical random walks lfl9l I20I1 . Here, we 
focus on those bipartite quantum walks defined in Refs. ifiol 
[l7ll for the purpose of obtaining quantum speed-ups in search 
problems. Such quantum walks, which we describe below, 
can also be derived from Ref. fTHl . 

To define the bipartite quantum walk, we first associate each 
classical state a of S with a quantum state \a) of an orthonor- 
mal basis of a d-dimensional Hilbert space ft. We then con- 
sider a tensor product Hilbert space Ha <8> of two copies 
of H. As in SA, we assume a given stochastic matrix M((3) 
describing the Markov process in S, with M((3) satisfying 
the detailed balance condition: TT a m a(7 i = 7r CT m a t a , with 
it a = e~@ E W/Z the components of the equilibrium distri- 
bution (Z = ^ e~P E ^ is the partition function). In the 
following we omit the dependence on (3 unless necessary. We 
define isometries X and Y that map states of H to states of 
Ha ® Hb as 

X\a) = \a)Y,V^W)i (2) 

a' 

Y\a') = ^2V^W)W) ■ (3) 

The symmetric operator H = X*Y, acting on H, has ele- 
ments h aa i = Jm^m^v fioll . Because of detailed bal- 
ance, we can write H = e^ H ^ 2 AIe~ ,3H ^ 2 , with H c the di- 
agonal operator H c \a) = E[a]\a). Therefore, the eigenvalues 
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A = 1 > Ai > ■ • • > X d -i > of H are those of M. If \<j>j) 
denotes the eigenstate of H with eigenvalue Xj, we have for 
J = 00] 



£i*> 



(4) 



The isometries X and Y define unitary operators Ux and 
Uy, acting on Ha ®Hb, via 



U x \ao) = X\a) 
U Y \oa) = Y\a) 



(5) 
(6) 



with |o) a selected state in H. The action of Ux and [7y in the 
remaining subspace is irrelevant. We now define Ri to be the 
reflection operator through the subspace spanned by {|(to)} 
and i?2 the reflection operator through the subspace spanned 
by {U x U Y \oa)}. Thus, 



i?i = 2IIi -1® 1, 

r 2 = 2n 2 - 1 ® i , 

where III and n 2 are the projectors 

n x = a® |o)(o| , 

n 2 = u x u Y (\o}(o\®i)u Y u x 



(7) 
(8) 



(9) 
(10) 



The unitary operation (rotation) W(M) = R2R1 defines the 
bipartite quantum walk based on the Markov chain M. This 
walk is related to the one used in Refs. Il0l[l7ll by a unitary, 
but /3-dependent, similarity transformation; using the trans- 
formed version is necessary for our QSA to work. 

The spectrum of W(M) can be directly related to the spec- 
trum of M ifioll . Defining the phases ip j = arccos A j , so that 



HI 



coscpjlfa) =X^Y\(f> j ) , 



(11) 



we have ipo = 0. When ip\ <C 1, the spectral gap of M (or 
if) is 1 - A x w {ipif/2. FromEqs. © and ©, 



H\U x U Y \o 4>j) = cosLp.j\4>j 0) 

IT2 1 <^ 0) = cos^j U x U Y \oi 



(12) 
(13) 



so the action of W(M ) in the (at most) two-dimensional sub- 
space spanned by 0) , U x U Y \o <fij)} is an overall Aipj ro- 
tation along an axis perpendicular to that subspace ||2lll . Thus 
the eigenphases of W{M) are ±2tpj, and its eigenvectors for 
j are: 



±i 



\/2 sin <pj 
When j = 0, we have 



(e^\^o)~U x U Y \o^)) 



\ip ) = \4>q 0) 



(14) 



(15) 



so a quantum algorithm that prepares the quantum Gibbs state 
I ?/)o) allows us to sample from the desired (equilibrium) dis- 
tribution by measuring Ha in the basis {|cr}}. All the other 



eigenphases of W(M) that were not described are either or 

7T. 

The (quantum) implementation complexity of Ux and Uy 
is proportional to the (classical) implementation complexity 
of a single step of the MCMC method described in Sec. [II] 
because U x , U x , U Y , and Uy may be implemented using 
a reversible version of the classical algorithm that computes a 
matrix element of M, It follows that the implementation com- 
plexity of W(M) is proportional to the classical complexity 
of implementing four steps in the MCMC method. 

The operations W(M) will be used below to implement 
the QSA. An important property that follows from our defini- 
tion of W(M) is that the overlap between the quantum Gibbs 
state \ipo(P)} and any other eigenstate in the 0-eigenphase 
subspace, at any (3 1 , is zero. To show this note that \<j>j) is 



a complete basis for H, and 
(j ^ 0). Thus, 



7I l 



d-1 



IMP)) = E*^» 



l^-i>] 



(16) 



3=0 



3=1 

Our algorithm uses this property to keep the state |^o(/3)) sep- 
arated form the remaining degenerate subspace. 



IV. QUANTUM SIMULATED ANNEALING ALGORITHM 

The QSA that we propose is basically a sequence of phase 
estimation algorithms (PEAs) projecting onto the quantum 
Gibbs state that is associated with the equilibrium state of S 
for different temperatures. The implementation complexity of 
SA is dominated by the gap of the stochastic matrix, which 
constrains the annealing schedule. For the QSA algorithm, 
the total implementation complexity is dominated by the im- 
plementation complexity of each PEA, given by the eigen- 
phase gap of the quantum walk. Because the latter is (quadrat- 
ically) larger than the former, the QSA algorithm results in a 
(quadratic) quantum speed-up of SA. 

We consider a sequence of inverse temperatures {(3k = 
kA/3}, with k = 1,...,Q, and p; = p Q = QAp. The 
choice of A/3 differs from the one used for SA. To understand 
the QSA, we begin by performing a Taylor series expansion 
of \MPk-i)) [Eq. ©] in p k . We obtain, 

\MPk-i)) = (1 - x ((jB>A ~ ^ c) ) lMPk)) 

+ 0{u 2 ), (17) 

where (E) 0k = J2 a E[a}e~^ E ^ / Z(p k ) = 
(4>o(Pk)\H c \(j)Q(Pk)) is the expectation value of the en- 
ergy (cost function), and v = A/3 Em- The (squared) overlap 
for two adjacent values of P is 



\(MPk)\MPk~i))\ 2 = 1 - o(v 2 ) . 



(18) 




|0>- 



S lo>- 
lo>- 



M io>- 

a io>- 



FIG. 2: Quantum simulated annealing algorithm. The algorithm is a 
sequence of Q calls to the PEA at /3i, . . . , j3f. After the last call, the 
state of TLa ® TCb is close to (/?/)) = l^o (/3/ ) o) , with proba- 
bility close to one. A measurement on Ha returns a state a in the 
ground state space of S with probability greater than 1 — e. 



FIG. 1: Phase estimation algorithm (subroutine) for the quantum 
simulated annealing algorithm. The first register of p qubits is used to 
encode the eigenphases of W(M (/3fe)). The second register denotes 
the bipartite system TLa ® T~Lb- The algorithm takes as input, in the 
second register, a quantum state sufficiently close to |^o(/3fc-i)}- A 
sequence of controlled W(M(/3k+i)) operations is performed and 
the inverse of the quantum Fourier transform is then applied; the 
composition of all these unitary operations is denoted PE(/3k). Fi- 
nally, the first register is measured. When the result of the measure- 
ment is such that the first register is projected onto |0) = |0i . . . P ), 
the PEA outputs a state close to \4>o{Pk)) in the second register. 



It follows that the probability of successful preparation of 
1 00 (/?/)}> a f ter Q = 0(\/v) projective measurements, can 
be bounded below by 1 — 0(y). This is called the quantum 
Zeno effect lfl4l[l5ll . Our QSA algorithm performs such pro- 
jections by calling the PEA at /3i, . . . , (3f. This technique was 
used in Ref. Wh to obtain the quadratic quantum speed-up for 
Graver's unstructured search problem. 

The PEA at the fcth step is depicted in Fig. [TJ The p ancil- 
lary qubits composing the first register are used to encode the 
eigenphases of W (M {(3k)) as binary fractions. In particular, 
2ipo = = [0i . . . 0p}2- The integer p is chosen to satisfy 
2 P = 0(l/(vv~5)). This choice allows us to bound the error 
due to the impossibility of representing the phases 2ipj with 
p bits (see Appendix lEl and Ref. lfl2ll ). The PEA gets as in- 
put a state close to |0 i[>o(f3k-i)}- It starts with a sequence 
of unitary gates that includes 2 P — 1 actions of the operation 
cW(M) = cR 2 cR l = u\U Y cP 0A U\U x cP 0B , controlled 
on the states of the first register (i = 1, . . . ,p). Here, 
cP 0A and cP 0B are the controlled selective sign change oper- 
ations on the states |o) of Ha and Hb, respectively. It contin- 
ues with an inverse quantum Fourier transform, and finally the 
first register is measured in the computational basis. For the 
given input state, the PEA outputs a state close to |0 ipoiPk)) 
with probability close to one. Since each use of cW(M ((3k)) 
has complexity proportional to that of four steps of the clas- 
sical MCMC method, the overall implementation complexity 
of the PEA is Npea = 0(l/(vVS)). 

The QSA is depicted in Fig. [2] It is composed of Q calls to 
the PEA, with a final measurement of Ha in the | cr) — basis. In 



Appendix|B]we show that, after the measurement, the proba- 
bility of finding Ha in the excited space can be bounded as 

V{a So) < de-Pri + r'Qv 1 , (19) 

for some constant r' = 0(1). We seek to make the above 
error of order e. Choosing (3f = 7 _1 log(2d/e 2 ), as in SA, 
makes the first term on the right hand side of Eq. ( fT9b of order 
0(e 2 ). Thus we need t'Qv 2 = 0(e). The condition QA(3 = 
(3 f implies A/3 = 0(e/((3 f E 2 M )) and Q = 0(((3 f E M ) 2 /e). 
Finally, because Mqsa = 0{QNpea), we obtain 

^^)=o((^) 3 M). 

(20) 

The above scaling with 1 /e 2 is for a single run of the QSA. 
Typically, repetition of the QSA makes the error exponentially 
low in the amount of resources used, so the dependence of 
Mqsa on e can be made logarithmic. The cubic scaling with 
the parameter Em/i is also worse than classical SA's linear 
scaling, but this is relatively unimportant as in most applica- 
tions this parameter will be bounded by a constant or a poly- 
nomial in instance size. 

Note that, since only the state of Ha is important for our 
purposes , the QSA can be implemented without measuring 
the ancillary qubits used in each PEA. In this case, the opera- 
tions FT -1 can be avoided 1I22I1 . This is because the quantum 
Zeno effect relies on the decoherence introduced by the inter- 
action with the ancillae, not the measurement itself. 



V. CONCLUSIONS 

We have presented a quantum algorithm to simulate classi- 
cal annealing processes by quantization of the simulated an- 
nealing algorithm implemented with MCMC methods. Such 
a quantization has been done by using techniques borrowed 
from quantum walks and quantum phase estimation. Our 
algorithm also exploits the quantum Zeno effect. We have 
shown that, if e denotes an upper bound to the probability of 
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not finding an optimal solution to a COP, the QSA requires 
resources Mqsa = O ^ log with ^ the spectral gap. 
Thus QSA outperforms SA in those problems where 5 <C 1, 
such as finding a ground state of a spin glass. SA requires 
■Ms A — 0(log(d/e 2 )/S) to assure the same error probability. 
Even if S A could be implemented more efficiently, the scaling 
of Ms a with <5 -1 may be unavoidable lfl8tl . Since initializing 
with a state close to \<f>o(f3f)) is not required by the QSA, our 
result has implications in the mixing time problem studied in 
Ref. HI. 

We expect that similar quantum speed-ups hold for the 
simulation of more general classical annealing processes. 
Moreover, our algorithm can easily be extended to simulate 
continuous-time annealing. Also, by choosing /3j = 1/T, 
with T > 0, the QSA can be used to speed up the calcula- 
tion of finite-temperature thermodynamic properties of classi- 
cal systems on a lattice. 

Finally, our QSA is one possible quantum algorithm to sim- 
ulate an annealing process. One may wonder if other quan- 
tum algorithms, based on quantum adiabatic evolutions, can 
still provide similar quantum speed-ups. The adiabatic theo- 
rem of quantum mechanics yields similar convergence rates. 
A simple, but not rigorous, proof is given by considering the 
adiabatic condition (cf. l24ll ): 



the uniform distribution). After P steps, this state evolves to 



d t /3(t) 



< < e, (21) 



with j 7^ 0. Other 0-eigenphase states have not been con- 
sidered as they do not overlap with \ipo(P)} at first order 
[Eq. (1161) 1. The overall implementation complexity of the 
adiabatic evolution (i.e., total evolution time) determined by 
Eq. ( f2Tb is 0(l/(e\/S)). Details will be given elsewhere. 
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APPENDIX A: CONVERGENCE OF CLASSICAL 
SIMULATED ANNEALING 

We now obtain an annealing schedule that assures conver- 
gence to the desired state when SA is implemented using 
discrete MCMC methods. The following analysis is based 
on Ref. [4], where similar rates have been obtained in the 
continuous-time case. Assume that we start with a state sam- 
pled from some probability vector fi(0) = 4(1, •■ • ,1) (i.e., 



0(0), 



(Al) 

with pk = kA/3. Because M is stochastic, normaliza- 
tion is preserved: Ecr=i H a (Pf) — 1. Let 7?(/3/) = 
(it 1 (f3 f ),..., ir d ([3 f)) be the desired (Boltzmann) equilib- 
rium distribution after the annealing process. That is, 

M(j3f)tf(j3f) = n(j3 f ), and also ELi^G 9 /) = L From 
the Cauchy-Schwarz inequality we obtain, for the probability 
of error, 



(A2) 



< 



\ 



E 



E 



Considering the worst case, in which all non-ground states 
have energy E[Sq] + 7 gives: 



E *°{p f ) < va t 



-Pn/2 



(A3) 



<J$So 



where 7 = min (7 g§ \E[cr] — £7[So] I is the spectral gap of E 
and d is the dimension of the state space S. Equation dA3b 
was obtained considering the worst case scenario in which the 
space of states having energy E[Sq) + 7 is highly degenerate. 
Thus 

?V (P) 2 §o) < Vde-^/ 2 \\h(P f )\\ 2 , (A4) 
where ||/i(/3/)||2 denotes the 2-normof 

_ ( ^(M /(/?/) \ 



h(j3 f ) 



(A5) 



To bound ||/i(/3/)||2, we define, as in Sec. Ill, the symmetric 
matrix H{(3 k ) = e^ H "/ 2 M(p k )e- l3kH ^ 2 , with H c the diag- 
onal matrix having E[l], . . . , E[d] as elements. We denote by 
AiOSfc) = 1 > A 2 (/3 fc ) > ■ ■ • > X d {Pk) > the eigenvalues 
of M(/3k) and H{f3k)- The eigenvector of H(Pk) with largest 
eigenvalue is 0] 



-d*( 



-P k E[l]/2 



(A6) 



where Z = E<x=i e~@ kE ^ is the partition function. Denote 
now as S — miiifc{l — \2(Pk)} the minimum spectral gap 
of the matrices H(f3f.) (° r ^I{Pk))- We will show that, when 
S <C 1, an annealing rate A/3 satisfying 



A/3E M < tS, 



(A7) 
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implies \\h(Pf)\\ 2 < V% JH- Here, E M = max CT \E[a}\ and 
r is a C(l) constant. 
We start by writing 

Afl(p k ) = ft{fi k+1 ) - fl(p k ) (A8) 
= (M(/3 k+1 )-l)fl(f3 k ) , 

where ^(/3 fe ) = nt=i M(/3 k ,)p,(0). Also, from the Taylor 
series expansion of n(/3 k ) and using Eq. (IA7t . we obtain 



V(«-vW) = 

= -H c )^W^) + 0{5 2 ) , 



(A9) 



where (JS)^ = T,t=i E[a]e-P kE W/Z is the expectation 
value of E at j3 k . Combining Eqs. ( IA8b and (IA91 . and defining 



we have 



Ah(p k ) = h((3 k+1 )~h((3 k ) 

= (H((3 k+1 )-l)h(f3 k )(l + (D(S)) 

- (±Ap({E) 0k - H c )^j h(p k ) + 0{5 2 ) . (A10) 

Therefore, if (• , •) refers to the standard inner product, 

(h((3 k ),Ah((3 k )) = 

= (h(p k ),(H(P k+1 )-l)h(P k ))(l + 0(S)) 



±A0{h(J3 k ), ((E) th - H c )h(h)) + C 



2 ) 
(All) 



The first term in Eq. (1A1 lb can be bounded by expanding 
h(ftk) as a sum of the eigenvectors of H(f3 k+1 ), denoted as 
{ej{(3 k+1 )}, with ei(/3fc+i) = 0r(/3 fc+ i) [see Eq. (|A6)]. 
Then, 

(H{p k+1 ) - l)h(/3 k ))(l + 0(S)) 

<-5(\\m)\\l-l)+0{5 2 ). 

(A12) 

This results in 

(h(/3 k ), Ah(/3 k )) 
< (-8 + ^Af3E M ^j \\h(0 k )\\ 2 2 + S + 0{5 2 ) , (A13) 

where we considered that (h(f3 k ) , H c h(f3 k )) < EM\\h(f3 k )\\ 2 
and, with no loss of generality, (E)/3 k > 0. Therefore, the 
increment on ||/i(/3fc)||2 can bounded as 

^h(ft k )\\l = \\h{p k+1 )\\l-\\h(j3 k )\\Z 

= 2(h(f3 k ),Ah(/3 k )) + \\Ah((3 k )\\ 2 2 (A14) 

< (-28 + A(3E M )\\h(Pk)\\l 

+ 2S + 0(5 2 ). (A15) 



Since ||/i(/3fc)||2 > 1 we have, for a proper choice of r = 0(1) 
in Eq. (TATt . 



|2 / x\\l.fa. MI2 , 2S 



A\\h(p k )U < -6\\h(p k )\\ 



Equivalently, 



\h((3 k+1 )\\ 2 2 <(l-6)\\h(p k )\\ 
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(A16) 



(A17) 



Furthermore, the condition 7?(0) = /2(0) yields to ||/i(0)||2 
1. Iterating Eq. (IA17b for k' = 0, . . . , k, we obtain 



\\h((3 k )\\l<2-(l-S) k <2 



(A18) 



Finally, using Eq. (IA4I) . we obtain the desired bound on the 
probability of error, given by 



(A19) 



APPENDIX B: IMPLEMENTATION COMPLEXITY OF THE 
QUANTUM SIMULATED ANNEALING ALGORITHM 

We first show how the PEA works for the eigenphases 

±2<fj of W(M), with ip = < <pi < ■ ■ ■ < (fd-i < tt/2. 
We write 

v 

2^ = 2n([.a{ . . . aj] 2 + Q) = 2tt(^ (Bl) 

i=l 

with \Q\ < 1/2 P+1 and 2n[.a{ . . . a^] 2 the bestp-bit approx- 
imation to 2ipj. The PEA (Fig. [TJ begins by applying a set 
of Hadamard gates to the p qubits in the first register, initial- 
ized in the state |0) = |0i . . . p ). These qubits are used to 
encode the eigenphases as binary fractions at the end of the 
PEA. The PEA then applies a set of operations W 2 ' (M), 
with i = 1, . . . , p, controlled on the states | 1,) of the first reg- 
ister. Consider the case where the initial state of TCa <8> 'Hb is 
one of the eigenstates \ip±j) of W{M) [Eqs. ( [Pfl ) and (fT5bl. 
The evolved joint state is 



^(|0 1 ) + e ± «°^)|l 1 ))--- 

■■•(|0 P ) + e ±ia '" 1 ( a ^|l p »|^) 



(B2) 



The next step is to apply the inverse of the quantum Fourier 
transform, denoted by FT^ 1 in Fig.Q] to the first register. Its 
action is given by 



2 P -1 



FT- l \m) = --= V 
\/9p — ' 



,-i2nmm'/2*\ I 



m'), (B3) 



where m, m' € [0, . . . , 2 P — 1] are natural numbers whose bi- 
nary representation denotes the states of qubits 1, . . . ,p. The 
evolved (joint) state is now 



^ 2 P -1 

to) = 2? 



-i2-nmm' /2 P ±im' (2^j ) 



m^ ±j ). (B4) 



m— 0,m'— 
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The final step of the PEA is to perform a projective 
measurement of the first register in the (computational) 
{|0j), basis (i = 1, . . . ,p). The probability of pro- 

jecting the first register onto some state |m) is determined by 
\o±j, m \ 2 , with 



°±3,: 



{mil>j\rj) 

2 P -1 
Op jL^i 



-i27rmm' j1 v Am' (2ifj ) 



\ \ _ e i[2 p (2ipj)-2mn] 



,,i(2ipj—2irm/2P) ' 



2P 1 



(B5) 



In particular, Oq m = 5 m and, since |1 — e lx \ > 2\x\/tt, we 
have |o-y m= o| < 7r /(2 p (2<^ J )). The error is due to the fact 
that, in general, 2ipj does not admit an exact representation 
using p bits. 

Clearly, the implementation complexity Afp ea of the PEA 
is of order 0(2 P ). The choice of p depends on the over- 
all probability of error of the QSA. Below we show that, by 
choosing |o±j >m= o| = 0(y), with v = A@Em, the QSA 
is guaranteed to succeed with a probability of error of order 
0(e). Furthermore, since mmj^{<pj(0)} — 0(y/S), where 8 
is the minimum spectral gap of M(/3), it is enough to choose 
p such that 2 P = 0(l/{v^fS)), giving a implementation com- 
plexity for each phase estimation Npea = 0(l/(vVS)). 

To obtain the implementation complexity of the QSA, it 
is helpful to consider the equivalent case where non of the 
measurements are actually performed until after the final 
PEA JH. The input state to the first PEA is \0 1 Vo(0)), where 
we introduce the subscripts l, . . . , q to denote the sets of p 
qubits used as ancillae in each PEA. The first PEA is per- 
formed at inverse temperature fi\ . From Eq. (JTVj 

10^0(0)) = (l-0(^ 2 ))|0 1 Vo(/3i))+OH|0 1 ^(/3i))- 



Also [Eq. ©J, 



(B6) 



d-l 



Wo-Wi)) = E -hl^+M) + CB7) 

After the implementation of the unitary PE{(3\) (see Fig. |2j, 
the evolved state is 

lo.VoOSi)) (B8) 

j,rn V 

Since only the states with m± — in the above sum contribute 
to the final probability of projecting onto |0i) at the end of the 



first PEA, it is convenient to rewrite Eq. ( IB8b as 

(i - o(v 2 )) lo, MPi)) + 0(^)10, 4>o(M) 

+ Q{v)\ X i). (B9) 

Here, (tp (f3 1 )\tp^(p 1 )) = (O^Xi) = and the order of the 
second term follows from the previous choice of p so that 

\o±j,m=o\ = 0(U). 

We now introduce the state |0 2 ) for the second set of p 
qubits, and evolve with the action of PE(^2). The output 
of the second phase estimation gives [Eq. (fTTI ll 

(i - eV)) |0z0x MPi)) + o(^)|o a o x ^(&)> 

+ 0{v 2 )PE(f3 2 )\QA ipo(Pi)) + 0(p)\ X 2) , (B10) 

with <0 1 a |x2> = 0. 

We repeat this procedure by introducing the states 
|0 3 ), . . . , |0„) and by evolving with PE(fo), PE{f3 Q = 
(if), respectively. Denote by |£) the evolved (joint) state of all 
the registers l, . . . , q and Ha After the measurement 

on l, . . . , q, the probability of projecting onto |0 q . . . Oi) is 
given by V Q = (£|P |£>, with P = |0 q . . . 0,) (0„ . . . OJ the 
projector onto the corresponding subspace. By a similar anal- 
ysis as the ones performed above for the first two steps, we 
obtain 



p \0^(i-o(v 2 )) Q \o,...o 1 MPf))+ 

O-i 

0{v 2 )P Q J2 PE{p Q ) ■ ■ ■ PE((3 Q ^ +1 )\0 q ... 0, ^(pQ-i)). 



(Bll) 



Thus the probability of Ha ®Hb being in the desired state 
\ipo(Pf)) can be bounded below, by using Eq. dB lib , as 



(l-0(v 2 )f -{Q-l)0(v 2 ) 



> 1 - t'Qv 2 



(B12) 



for some constant t' = 0(1). 

Assume now that the state of Ha ® Hb is \4>o((3f)) = 
\(po{Pf)o) = J2t=i V ' 7r ° 'f)\° °) ■ If a measurement on the 
| a)— basis is performed on Ha, the probability of finding the 
system in an excited state can be bounded by de~° 1 "' ' , with 7 
the spectral gap of E. Thus, after the QSA, the total probabil- 
ity of such an event, which is the error probability for QSA, 
can be bounded above by 



V{<7 & S ) < de 



-Pn 



(B13) 



as claimed. 
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